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Abstract 

The ocean is a complex, constantly changing, highly dynamical system. Prediction 
capabilities are constantly being improved in order to better understand and forecast 
ocean properties for applications in science, indnstry, and maritime interests. Our 
overarching goal is to better predict the ocean environment in regions of complex 
topography with a continental shelf, shelfbreak, canyons and steep slopes using the 
MIT Multidisciplinary Simulation, Estimation and Assimilation Systems (MSEAS) 
primitive-equation ocean model. We did this by focusing on the complex region 
surrounding Taiwan, and the period of time immediately following the passage of 
Typhoon Morakot. This area and period were studied extensively as part of the 
intense observation period during August - September 2009 of the joint U.S. - Taiwan 
program Quantifying, Predicting, and Exploiting Uncertainty Department Research 
Initiative (QPE DRI). Typhoon Morakot brought an unprecedented amount of rainfall 
within a very short time period and in this research, we model and study the effects 
of this rainfall on Taiwan’s coastal oceans as a result of river discharge. We do this 
through the use of a river discharge model and a bulk river-ocean mixing model. We 
complete a sensitivity study of the primitive-equation ocean model simulations to the 
different parameters of these models. By varying the shape, size, and depth of the 
bnlk mixing model footprint, and examining the resnlting impacts on ocean salinity 
forecasts, we are able to determine an optimal combination of salinity relaxation 
factors for highest accuracy. 

Thesis Supervisor: Dr. Pierre F. J. Lermusiaux 
Title: Associate Professor of Mechanical Engineering 
Massachusetts Institute of Technology 
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Chapter 1 


Introduction 


“The seas unite the world, no matter if individual nations are landlocked or not. Our 
individual sovereignties, our freedoms, and our prosperity are linked by the eommon 
medium of the oceans. ” 

-Admiral Mike Mullen 


The ocean has an enormous impact on our lives, and yet relatively little is known 
about it. The ocean is a complex, constantly changing, highly dynamical system. 
Harnessing the knowledge of the depths gives the advantage to the owner of that 
knowledge, whether his purpose is scientific research, navigation, fishing, acoustic 
communications, shipping, climate change, or naval warfighting. 

In this work, our goal is to better predict the ocean environment in a region of 
region of complex topography with a continental shelf, shelfbreak, canyons and steep 
slopes. We add to the challenge the ocean’s response to an extreme weather event. 
Today’s ocean models resolve complex physical processes taking place on a wide range 
of spatial and temporal scales. Understanding the coupled ocean-atmosphere system 
and enhacing the performace of our models allows a reduction in uncertainty- an 
important factor in making the right decisions at the right time. The ability to 
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reduce uncertainty in our knowledge of the surrounding environment, in both present 
and future situations, is the ability to gain the advantage. 

Specihcally, we study the region around Taiwan, in the aftermath of Typhoon 
Morakot in August 2009. The Quantifying, Predicting, and Exploiting Uncertainty 
Department Research Initiative (QPE DRI) seeks to answer these questions and more 
with a targeted observation period and modeling studies. Part of this DRI involved 
the coupling of observations with modeling studies, using the sophisticated ocean 
models from the Multidisciplinary Simulation, Estimation and Assimilation Systems 
(MSEAS) group at MIT. In the present research, we study the effects of the unprece¬ 
dented rainfall on Taiwan’s coastal oceans as a result of river discharge. We do this 
through the use of a river discharge model and a bulk river-ocean mixing model, and 
we rehne that model to find the combination of input parameters that most closely 
represent observed ocean properties. Specifically, we complete a sensitivity study 
of the primitive-equation ocean model simulations to variations of the bulk mixing 
model parameters in the coastal regions directly adjacent to the river mouths. 

Chapter two discusses the background information used for this research, a sum¬ 
mary of the research initiative utilized in the project, and a description of the river 
discharge and river-ocean mixing models that we seek to improve upon. This chapter 
also includes a review of some of the literature that is applicable to this research. 
Chapter three describes the methods used and further developed in this research, 
outlining some of the steps taken and the reasoning for them. Chapter four describes 
the results obtained in the course of this research. Chapter five concludes the work 
and details possible future work. 
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Chapter 2 


Background and Literature Review 


This chapter discusses the Taiwan region and complex oceanography, then provides 
an explanation of past projects used to study the region. The current state of the 
river discharge and bulk mixing model is examined. 


2.1 Taiwan 


Taiwan is an island in the northern west Pacific Ocean, about 180 km (100 miles) from 
the southeast coast of China. It is bordered by the Philippine Sea to the east, the 
Luzon Strait and South China Sea to the south, the Taiwan Strait to the west and the 
East China Sea to the north. The total landmass is 35,980 sq km. The eastern part 
of the island is dominated by mountains, while the western part is flatter, and home 
to the majority of the island’s 23 million people (U. S. Central Intelligence Agency, 
2014). The climate is tropical, dominated by the winter northeasterly monsoons and 
the weaker summer southwesterly monsoons (Rudnick et ah, 2011). On average, 3 
- 4 typhoons make landfall each season. The island lies along an active fault, and 
experiences frequent earthquakes (Taiwan Central Weather Bureau, 2014). 
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Figure 2-1: Bathymetric chart and topography in the seas surronnding Taiwan. (Rud- 
nick et ah, 2011) 


2.1.1 Taiwan Oceans and Topography 


The oceans surronnding Taiwan have a highly complex underwater topography. The 
island was formed when the Philippine Sea Plate subducted under the Eurasia Plate, 
(Jan et ah, 2011) along the Ryukyu Trench (Rudnick et ah, 2011). To the west of 
Taiwan is the continental shelf, and to the east is a steep continental slope, ending 
in the Huatnng Basin, 6,200 m deep at its lowest point (Rndnick et ah, 2011). A 
series of canyons (the Keelung Valley, the Mien-Hua Canyon, and the North Mien- 
Hua Canyon) lie to the north of this basin, cutting across the shelf (Jan et ah, 2011). 
To the east is a series of trenches and ridges, inclnding the 1-Lon Ridge (Jan et ah, 
2011). South of the Huatung, the Luzon Island Arc is the top of a series of nnderwater 
mountains, extending south to the larger Philippine Islands (Rudnick et ah, 2011). 
See Fignre 2-1. 
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2.1.2 Taiwan Ocean Features 


The region’s oceanographic features are dominated by the Kuroshio current offshore 
to the east; an area of intermittent upwelling over the continental slope commonly 
referred to as the "cold dome" forms off the northeastern tip of the island (Tang et ah, 
1999). 

As presented by Rudnick et ah (2011), the Kuroshio current is a western boundary 
current of the north Pacific subtropical gyre (see Figure 2-2). The Kuroshio, along 
with the westerly North Equatorial Current (NEC) and the southerly Mindanao cur¬ 
rent, forms the major circulation pattern in the northwest Pacific. The Kuroshio and 
the Mindanao split from the western edge of the NEC at the Philippine Islands as 
shown in Eigure 2-3. The water masses making up the NEC are the North Pacific 
Tropical Water (NPTW), located by a salinity maximum at 200 m, and the North 
Pacific Intermediate Water (NPIW), with a salinity minimum at 500 m. Strong eddies 
form off the western boundary, near the coast of Taiwan. 

The Kuroshio, like other north-south currents, is relatively narrow and swift, 
ranging between 80-100 km, and having maximum speeds between 1.5-5 knots (Talley 
et ah, 2011). The current, measured just northeast of Taiwan, has a mean transport 
of 21 Sv (Rudnick et ah, 2011). 

In the region surrounding Taiwan, specific seasonal changes in the current’s path 
occur, shown in Figure 2-4. During winter, the Kuroshio is affected by the seasonal 
northeasterly monsoon and some Kuroshio water fiows over the continental shelf. In 
summer, the monsoon shifts and the Kuroshio moves northeastward off the continental 
shelf. Local oceanographic properties are substantially impacted by this seasonal shift 
(Jan et ah, 2011). 

During the summer months, the Kuroshio is no longer opposed by the winter 
monsoon and turns away from the northeast corner of Taiwan. This allows the ocean 
north of Taiwan to be dominated by Taiwan Strait waters, which are fresher and 
warmer than the surrounding waters. The passage of typhoons can disrupt this 
pattern, allowing cold, salty Kuroshio water to be brought to the surface via upwelling 
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2-2: Surface currents of the Pacific Ocean, shown during winter monsoon 
(Talley et ah, 2011) 
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Figure 2-3: Origins of the Kuroshio current. (Rndnick et ah, 2011) 
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Figure 2-4: Seasonal variations in the Kuroshio near Taiwan. (Rudnick et ah, 2011) 


over the shelf. This water mixes with the shelf water, creating a relatively short¬ 
lived effect known as the “cold dome.” During the winter, the cold dome structure 
is masked by intrusions from the Kuroshio and the lack of horizontal temperature 
variations caused by seasonal atmospheric cooling. (Jan et ah, 2011) 

We now summarize the work of Jan et al. (2011). The cold dome is centered at 
25.625°N, 122.125°E, with a diameter of 100 km. The average temperature at the 
center is 21°C. In the region, the temperature standard deviation at 50 m is 2.5°C, 
while in the cold dome the standard deviation is only 2°C. At 100 m, the cold dome 
has a standard deviation of <1°C. These small standard deviations indicate that the 
cold dome waters are sourced by upwelling from below 100 m. ADCP data indicate 
a weak cyclonic flow at 50 m, another sign of upwelling within the dome structure. 

The oceanographic features of the region are impacted by many forcings, includ¬ 
ing atmospheric forcings (especially typhoons and monsoons), buoyancy forcing from 
freshwater river discharge, eddies, and water mass formation. The presence of smaller 
currents, wave propagation, and internal tides all contribute to a highly complex ocean 
environment that is challenging to predict (Gawarkiewicz et ah, 2011). 
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Figure 2-5: Three-dimensional rendering of cold dome formation. (Jan et ah, 2011) 

2.2 Quantifying, Predicting, and Exploiting Uncer¬ 
tainty Department Research Initiative 

The Quantifying, Predicting, and Exploiting Uncertainty Department Research Ini¬ 
tiative (QPE DRI) was undertaken in order to study the “wide variety of physical 
processes ocurring on a broad range of spatial and temporal scales'^ in the region 
snrrounding Taiwan (Gawarkiewicz et ah, 2011). 

2.2.1 Motivation for and Goals of the QPE DRI 

Model forecasts are crucial to the exploration and understanding of the ocean. How¬ 
ever, model forecasts will inherently have some uncertainty, and being able to measnre 
nncertainty is important to the users of model forecasts. This is especially important 
in a complex environment such as the ocean surrounding Taiwan. 

The QPE DRI was a joint US-Taiwan program aimed at exploring uncertainty 
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and forecasting in the highly dynamical and challenging environment around Taiwan. 
Led by Glen Gawarkiewicz of the Woods Hole Oceanographic Institution and Jan Sen 
of the National Taiwan University, the program included participants from multiple 
organizations within the US and Taiwan, and was partially sponsored by the OlRce 
of Naval Research. 

One of the goals of the QPE project was to better understand the complex dy¬ 
namics of the region by combining model forecasts with in-situ observations and by 
performing sensitivity studies on existing models. By collecting data and carefully 
processing it, we can gain an increased understanding of the ocean processes, as well 
as improved model forecasts. Uncertainty within model forecasts must be quantified 
in order for the model forecasts to be useful. 

Other program goals, though not neccessarily discussed in this thesis, are to couple 
acoustic and environmental models with tides and nesting, research more effective 
ways to link small-scale and large-scale models, and use adaptive sampling techniques 
(Lermusiaux et ah, 2010). 

The late summer was chosen as the timeframe because the unique combination 
of the summer monsoon, typhoon frequency and cold dome formation provide ample 
opportunities to study the dynamic oceanography within the region (Gawarkiewicz 
et ah, 2011). 

2.2.2 lOP Observation Plan 

The Intense Observation Period lasted from 18 August to 10 Septemeber 2009. Gliders 
were released into the Kuroshio in May 2009 by the US research vessel R/V Roger 
Reville, and near surface drifters were released weekly in the leadup to the lOP. This 
method obtained observations for model initialization and later comparison. 

From 13 - 17 August, eighty-eight Surface Velocity Program (SVP) drifters were 
deployed; the Taiwanese research vessels R/V Ocean Researcher 2 (OR2) and R/V 
Ocean Researcher 3 (OR3) launched half of these into the waters northeast of Taiwan 
in early August, and the R/V Roger Reville launched the other half during the lOP 
(Gawarkiewicz et ah, 2011). 
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Figure 2-6: Map of observing plan for QPE lOP. Triangles indicate CTD drops of R/V 
Ocean Researcher 2 (OR2) and squares indicate CTD drops of R/V Ocean Researcher 
3 (OR3). Yellow star is a biogeocliemcial sampling station, and gray arrows show 
mean current velocity, taken from historical ADCP data. Red dashed line shows the 
track of Typhoon Morakot. (Jan et ah, 2013) 

The OR2 and OR3 also took measurments via Conductivity-Temperature-Depth 
sensors (CTD) and Acoustic Doppler Current Profilers (ADCP) during this period, 
and again from 21 - 27 August and 27 August - 2 September (Jan et ah, 2013). The 
locations for these profiles are shown in Figure 2-6. 

Two SeaSoar surveys near mooring arrays to the east of Taiwan and profiler fioats 
over the continental shelf provided high-resolution measurements. SeaSoar is a towed 
vehicle that can undulate, allowing an oceanographic profile to be measured approx¬ 
imately every 3 km (Woods Hole Oceanographic Institution, 2014). 

Bottom-mounted ADCPs and CTDs were also used in the continental shelf region 
and the North Mien-Hua Canyon (labeled NMHC in Figure 2-6) (Gawarkiewicz et ah, 
2011 ). 

During the lOP, the Multidisciplinary Simulation, Estimation and Assimilation 
Systems (MSEAS) model simulations were run 10 - 30 times daily, incorporating new 
data into the model almost as soon as it was collected (Haley et ah, 2014). The 
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real-time updates allowed adaptive sampling plans to be utilized for both physical 
and acoustic sampling. This meant the sampling plan could be guided by the model 
forecasts, and allowed the model forecast uncertainties to be improved by utilizing 
the best possible sampling plan. 

2.3 MSEAS Primitive-Equation Ocean Model 

The goal of the MSEAS primitive-equation ocean model is to produce accurate and 
timely regional ocean forecasts, accounting for tidal-to-mesoscale oceanic processes. 
The model specializes in regions with complex and varied topography. It is capable 
of implicit two-way nesting, resolving multiscale hydrostatic primitive equations with 
a nonlinear free surface, or high-order hnite element code on unstructured grids for 
non-hydrostatic processes (Haley et ah, 2014). 

The MSEAS family of models includes initialization schemes, nested data-assimilative 
tidal prediction and inversion; fast-marching coastal objective analysis; stochastic 
subgrid-scale models; generalized adaptable biogeochemical modeling system; La- 
grangian Coherent Structures; non-Gaussian data assimilation and adaptive sampling; 
dynamically-orthogonal equations for uncertainty predictions; and machine learning 
of model formulations (Haley et ah, 2014). 

We summarize, from Haley and Lermusiaux (2010), the primitive equations used 
in the MSEAS ocean model. These are derived from the Navier-Stokes equations and 
assume hydrostatic and Boussinesq approximations. 


Cons. Mass 

_ ^ dw 

V • u + = 0 , 

oz 

(2.1) 

Cons. Horiz. Mom. 

Du pt 1 „ 

— + fkxu = - Vp + F , 

Dt po 

(2.2) 

Cons. Vert. Mom. 

dp 

(2.3) 

Cons. Heat 

DT _ T 

Dt ^ ’ 

(2.4) 

Cons. Salt 

pS 

Dt ’ 

(2.5) 

Eq. of State 

P = p{z,T, S), 

(2.6) 
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Eq. FreeSurface 


0 


(2.7) 



As noted in Haley and Lermusiaux (2010), the state variables are the horizontal 
and vertical components of velocity {u,w), the temperature, T, and the salinity S. 
The ^ is the 3D material derivative, p is the pressure, / is the Coriolis parameter, p 
is the density, po is the (constant) density from a reference state, g is the acceleration 
due to gravity and k is the unit direction vector in the vertical direction. The gradient 
operators, V, in Equations 2.1 and 2.2 are two dimensional (horizontal) operators. 
The turbulent sub-gridscale processes are represented by F, and . Finally, the 
surface elevation is g and H = H{x,y) is the local water depth in the undisturbed 
ocean. 

For the equations governing the numerical methods, as well as the numerical 
discretizations utilized in the model, see Haley and Lermusiaux (2010, Sect. 2.1 and 
2 . 2 ). 

2.4 Typhoon Morakot 

Typhoon Morakot was a Category I typhoon that affected Taiwan and the surrounding 
region in early August, 2009. 

2.4.1 Typhoon Morakot Overview 

Typhoon Morakot developed over the Philippine Sea on 04 August 2009 and intensi¬ 
fied into a typhoon by the next day as it tracked westward. Morakot had maximum 
sustained winds of 80 kts, making it a Category I typhoon (Joint Typhoon Warning 
Center, 2009) as it approached Taiwan at a speed of about 11 kts (20 km/hr) (Jan 
et ah, 2013). The storm spread for a diameter of more than 1,000 miles (1700 km) 
(Gutro, 2009). The size of the storm is shown in Figure 2-7. The storm eye moved 
over Taiwan just before 1800Z on 07 August 2009 (Joint Typhoon Warning Center, 
2009) and interaction with Taiwan’s rugged topography caused its speed to decrease 
by half (Jan et ah, 2013). The storm eye moved offshore towards China before 0600Z 
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Figure 2-7: Satellite image of Typhoon Morakot on 06 August 2009 at 0525Z, taken by 
Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the NASA Aqua 
satellite. (NASA/MODIS Rapid Response, 2009) 


on 08 August. 


2.4.2 Precipitation from Typhoon Morakot 

The storm center passed over the island in less than twelve hours; however, the slow 
speed, immense span, and intake of moisture combined for the highest rainfall totals 
in Taiwan in 50 years. NASA’s Tropical Rainfall Measuring Mission satellites showed 
nearly half the island received more than 600 mm during 03 - 10 August 2009, and 
some areas showed rainfall totals of greater than 1000 mm (40 inches) (Lang, 2009). 
These estimates are lower than those of the Taiwan Central Geological Survey, which 
had high-precision instruments with in-situ measurements throughout the period. 
Their measurements show a large portion of southern Taiwan received in excess of 
2.0 m of rainfall, with some areas receiving 2.6 m over the period of the storm. See 
Figure 2-8. 
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Total COAMPS Precipitation (mm) - 03-10 August 2009 



(a) (b) 



(c) 


Figure 2-8: Precipitation estimates for rainfall associated with Typhoon Morakot ob¬ 
tained from (a) NASA’s Tropical Rainfall Measuring Mission satellites (SSAI/NASA, 
Hal Pierce, 2009), (b) nowcast from the NRL’s Coupled Ocean/Atmosphere Mesoscale 
Prediction System (COAMPS) (Naval Research Laboratory, 2012), and (c) the Tai¬ 
wan Central Geologial Survey in-situ measurements (Liu, 2010). 
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The abnormally high rainfall contributed to mudslides and severe flooding (Lang, 
2009). Taiwan’s rivers, unable to absorb the excess precipitation, showed record 
discharge rates, peaking 6-10 August. 

The typhoon had devastating consequences to the island of Taiwan, as well as 
mainland China. For oceanographic research, the timing of the storm allowed for a 
rare opportunity to compare oceanographic data before and after the passage of a 
typhoon. 

2.5 River Discharge Model 

The unprecedented amount of precipitation was unable to be absorbed by the ter¬ 
rain. Due to lack of time and personnel during the DRI, the river runoff from the 
ffooded rivers was not incorporated into the real-time MSEAS ocean model. However, 
the freshwater river discharge must be accounted for in ocean models to accurately 
forecast oceanographic variables, including salinity. Ocean models use forecasts from 
atmopsheric models, but in this case, precipitation was underforecast, so the direct 
impact on the ocean’s properties could not be accurately represented. Nonetheless, 
the majority of the freshwater input was later found to come from the river dis¬ 
charges, in part because of the higher level clouds encountering the mountain range 
(see Section 2.1). To properly input the river component of this additional freshwater 
source, a river discharge model was created (Mirabito et ah, 2012). The following is 
a summary of the development of the river discharge model; for the more detailed 
explanation of the model, we refer to the original description by Mirabito et al. (2012). 

2.5.1 River Discharge Model Concept 

A control volume is defined such that precipitation enters the grid box through the top 
boundary. Some loss due to soil absorption and collection in reservoirs occurs at the 
bottom boundary. All other precipitation leaves via lateral boundaries as freshwater 
river discharge (see Figure 2-9). Defining the control volume as fl with a boundary 
dfl, and using the traditional notation of p for rainwater density in kg/m^ and u for 
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Figure 2-9: Schematic of river discharge model control volume. 


flow velocity in m/s, it is shown that: 


[ ^ dQ + [ pu ■ fi dA = 0. 
Jn ot Jdci 


( 2 . 8 ) 


Assuming incompressible rainwater with constant temperature and salinity allows 
the first term in Equation 2.8 to become 0. The second term is expanded to each side 
of the control volume. 



pu • fi dA = 


~Qp ++ Os — 0) 

i 


(2.9) 


where Qp is the precipitation rate (m^/s), Qi is the total discharge rate (m^/s) at 
river mouth i, and Qg is the total seepage rate (m^/s). Evaporation is neglected 
in Equation 2.9. The river discharge model assumes that some of the rainfall will 
be absorbed into the soil or collect in lakes and reservoirs. For simplicity, this is 
assumed to be a constant fraction, neglecting changes in terrain and soil properties. 
Per estimates from the National Taiwan Normal University, one-third of the rainfall 
is assumed to be lost to soil seepage, irrigation or collection in natural lakes and 
reservoirs, such that Qs = aQp where a = |. Using Qi to represent the total river 
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discharge for one river i during the period: 


^4 = (i-«)0p. (2.10) 

i 

To apply this concept to the MSEAS ocean model, the freshwater source is mod¬ 
eled as a salinity sink. A low-salinity input is applied to selected gridpoints near 
the river mouth. The execution is more complicated, due to the complexity of the 
involved processes and the lack of data. 

2.5.2 Precipitation Estimation and Drainage Basins 

The data for rainfall during the period of heaviest rainfall, 6-10 August 2009, 
provided in Figure 2-8 provides total precipitation amounts, but not precipitation 
rates. The island was divided into four drainage basins, with boundaries based on 
information from the Water Resources Management Research Center at the National 
Taiwan Ocean University. The total precipitation for each basin, denoted Qpj, was 
estimated from Figure 2-8(c). Applying Equation 2.10 provides a total discharge 
amount per basin over the period: 


Q, = {l-a)Qpj (2.11) 

where the index j indicates the drainage basin. 

Now that Mirabito et ah (2012) have estimated for the discharge per basin, the 
next step is to compute the discharge associated with each of the major rivers. Taiwan 
has 129 rivers, and twenty-one major rivers. Only eleven had readily available, reliable 
annual mean discharge data; these rivers account for 53% of the total river discharge 
of Taiwan, and will be used in the model to represent 100% of the drainage from 
Morakot. See Figure 2-10. These eleven rivers were each assigned to a drainage basin 
based on their geography. Using the annual mean discharge rate for each river, a 
weighting factor /5j can be computed for each river i, which gives an approximation 
of the distribution of discharge within the basin. This allows the total discharge 
amount per basin Qj to be divided appropriately among the rivers within Basin j. 
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(a) 


(b) 


Figure 2-10: Figures used in estimation of individual river discharge, (a) Map of major 
rivers of Taiwan, (b) Watersheds in Taiwan and model discharge basin boundaries. 
Numbers indicate mouths of major rivers, while shading indicates River Discharge 
Model drainage basins. (Mirabito et ah, 2012) 

In the following equation, rij represents the total number of rivers assigned to basin 

j- 



( 2 . 12 ) 


2.5.3 Time Dependency of River Discharge 

As discussed above, the data lends itself to computations of discharge amounts, but 
it would be physically unrealistic to input the full discharge amount into the ocean 
model over a single time instant. Rather, a discharge rate must be computed in order 
to realistically simulate the addition of the freshwater discharge into the ocean. Thus, 
an average discharge rate, denoted Qj, can be computed for each basin by taking the 
integral over the time period, with to representing the start of our data at OOOOZ on 
6 August and T representing the 96 hours over which the heaviest rainfall occurred: 
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Basin 
Index (j) 

River 
Index (i) 

River Name 
Discharge (m^/s) 

Annual Mean 

Pi 

1 

1 

Gaopi'ng (i^M) 

220 

0.81 


2 

Zengwen 

50 

0.19 

2 

3 

Zhuoshui (illzK) 

210 

0.51 


4 

Wu (B) 

120 

0.29 


5 

Dajia (^¥) 

50 

0.12 


6 

Daan (7^^) 

35 

0.08 

3 

7 

Danshui (ti^zK) 

200 

0.80 


8 

Lanyang (Ml^§) 

50 

0.20 

4 

9 

Hualian (iftS) 

60 

0.31 


10 

Xiuguluan (^M'M) 

55 

0.28 


11 

Bmnan (.^1^) 

80 

0.41 


Table 2.1: The eleven rivers of Taiwan used in the river discharge model and their 
assigned drainage basin. The relative strength of the annual mean discharge of each 
river within its basin determines the discharge weight /3j. (Mirabito et ah, 2012) 


Qi = 13000 m^/s Qa = 3500 m^/s 
Q 2 = 10000 m^/s Q 4 = 6300 m^/s 


Table 2.2: Average Discharge Rate by Basin. (Mirabito et ah, 2012) 


— Aof 1 /‘^0+r _ 

Q^ = - Qidt, zg{1,2,...,11}. (2.13) 

j- Jto 

This results in an average discharge rate per river for the period of heaviest rainfall 
during Typhoon Morakot. While more reasonable than inputting the entire discharge 
amount at once, this rate is still not the best physical representation, as the actual 
discharge rate will fluctuate highly during this time period. 

In order to best represent the freshwater discharge rate throughout the typhoon, 
Mirabito et al. (2012) uses a time-dependent scaling factor Aj which must be computed 
for each river i such that the best representation of the discharge rate per river Qi 
can be modeled: 


Qi — \Qref i G {1; 2,..., 11}. 


(2.14) 
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In-situ measurements of river discharge were available throughout the period for 
only two of the above rivers; the Zhuoshui River (index i = 3) had readings taken ap¬ 
proximately every 7.5 hours, and the Gaopmg River (index i = 1) had daily readings. 
These are, in terms of annual mean discharge, the largest two rivers on Taiwan. The 
data covered the period from 1 August to 15 September 2009. Mirabito et al. (2012) 
used this data to represent the average normalized discharge rate for each river. 

The known discharge rate data, denoted Qz and Qg, respectively, was extrapo¬ 
lated so each had the same timescale. It was then normalized to account for differences 
in discharge quantity, and the average of the two series was taken, denoted \Qref\- 
This average can then be used to represent the time series for discharge rates for the 
other nine rivers. 


2 {vg^Vz)' 

where 

Vg = f QGdt (2-16) 

Jti 

with representing 0600Z 1 August and T* representing 45 days, the length of the 
data set for Qz and Qg- 

It can be shown that the period of heaviest rainfall, lasting 96 hours, accounts for 
67.34% of the total rainfall during the 45-day period represented by Qref'- 



0.6734 


which allows computation of the scalar quantity A* via 


(2.17) 


rto+T 0.6734 /■‘o+i'a , 

/ Qref dt — rpn / Qi dt^ 
Jto 1 PiQi Jto 


(2.18) 


TfiiQj 

0.6734' 


(2.19) 


Now the average discharge rate Qi can be multiplied by the normalized time series 
to determine how much freshwater discharge to apply to each timestep in the model. 
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Peinan River 


Hualien River 


Hsiukuluan River 


Kaoping River 



Day (August 2009) 
Tsengwen River 



Day (August 2009) 
Taan River 




Day (August 2009) 
Choshui River 


Day (August 2009) 
Tanshui River 





Day (August 2009) 
Wu River 


Day (August 2009) 
Lanyang River 





Day (August 2009) 
Tachia River 



Day (August 2009) 


Day (August 2009) 


Day (August 2009) 


Day (August 2009) 


Figure 2-11: Discharge rates for each river as computed by this method, with the 
actual data replacing the time series for the Zhuoshui and Gaoping Rivers. (Mirabito 
et ah, 2012) 


The actual time series replaced the averaged time series for the Zhuoshui and Gao¬ 
ping Rivers to increase accuracy. Figure 2-11 shows the discharge rates for each river 
as computed by this method, with the actual data replacing the time series for the 
Zhuoshui and Gaoping Rivers. 

The river discharge model of Mirabito et ah (2012) is now able to provide a salinity 
sink of appropriate size and time series to model the freshwater discharge. 


2.5.4 Forcing the MSEAS Ocean Model with the River Dis¬ 
charge Model 

In the original work of Mirabito et al. (2012), the modeled river discharge was applied 
to the ocean model simulations using a salinity relaxation factor in the hnite-volume 
grid boxes nearest to the respective river mouth. For the initial tests, the size of the 
grid footprint was limited to seven grid cells and the salinity relaxation was applied to 
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a depth of 10 m. The river discharge computed by the above methods was distributed 
evenly throughout these seven gridpoints and throughout the depth. 

2.6 Preliminary Results of the Original River Dis¬ 
charge Model 

When applied to the MSEAS ocean model, near-shore salinity was signihcantly im¬ 
proved during and after typhoon passage, more closely matching the data collected 
during the QPE experiment. In Figure 2-12, the impact of the river discharge model 
on forecasted salinity is obvious. Several areas of improvement were identified and 
will be addressed in Section 3.2. 
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(d) August 17 - with river 
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(e) August 19 


(f) August 25 - with river 


Figure 2-12: Surface salinities forecast by the MSEAS ocean model with no river 
forcing (left column) and with river discharge model (right column). (Mirabito et ah, 
2012) 
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Chapter 3 


Methods 


In this chapter, we discuss the different approaches to improving the river discharge 
model and the bulk mixing model. The chapter begins with a discussion and deriva¬ 
tion of the salinity forcing equation, then moves on to specific areas of the models 
that were explored for possible improvement. These improvements have the benefit 
of hindsight, but incorporating these changes into future use will enable enhanced 
forecasting. 


3.1 Bulk Mixing Model of River and Ocean Waters 

The river discharge model and freshwater fiow-rate estimates provided in Section 2.5 
need to be applied as a forcing to the ocean model. The subsequent mixing of river 
and ocean waters is discussed in this section where we derive a simple mixing model. 

Normally, the river water enter the ocean in river mouths and estuaries, possibly 
with hydraulic jumps and a river head differential, usually forming small-scale plumes, 
vertical billows with instabilities such as Kelvin-Helmoltz instabilities and overall 
turbulent structures in these river mouths and estuaries (Ozgokmen and Chassignet, 
2002). Subsequently, complex advection pattern and additional eddy and turbulent 
mixing occur. 

In this section, we will represent all of these processes in the vicinity of the river 
mouths as a bulk volume mixing of input fresh riverwater with the local ocean water, 
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provided at a rate proportional to the local river discharge considered. This local 
mixing is modeled as a source, or forcing, of lower salinity river water that is input 
to the selected local ocean volume (specified by an ocean area and depth range near 
the river mouth). That given local ocean volume is itself discretized by the set of 
finite-volumes in the MSEAS ocean model, and the freshwater river water is provided 
in each of these finite-volume cells, usually at a rate equally distributed to each of 
these finite-volume cells. 

Next, we provide the derivation of this lower salinity source or forcing function. 
Then we look at possible areas of improvement in the original model, then explore 
these potential improvements. 


3.1.1 Derivation of Salinity Forcing Function 

Here we consider a local ocean finite-volume V selected near the mouth of a given 
river R. This volume, affected by the fresh water river forcing, is specified by an ocean 
area and depth range near the coast, in accordance with the expected ocean currents 
in the region. We assume that a source, or forcing, of lower salinity river water is 
input to this selected local ocean volume. In Figure 3-1, we sketch such a volume, 
both without and with the fresh water river discharge. In our notation, we define the 
mass flow rate of salt to be rfis = S -Q where S is the salinity concentration in kg/w? 
and Q is the flow rate of the water in irfi/s. 

We begin by applying the conservation of mass of salt to the considered finite- 
volume V. We obtain, using conservation in the summation form, 

ms,in - Y '^S,out- (3.1) 

We then rewrite the equation in terms of the definitions above to be 

j^Ms = Y (SQU - Y iSQ)out. (3.2) 

Consider a control volume, or grid cell, with no river discharge model. In this case, 
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the only salt fluxes are the flux into and the flux out of the cell via normal ocean 
processes. 

= SinQin — SoutQout (3-3) 

If we now consider the conservation of mass of water, assuming incompressibility at 
first order, i.e. div u = 0 in a Boussinesq fluid (Cushman-Roisin and Beckers, 2011), 
and neglecting the volume variations due to the free surface (Haley and Lermusiaux, 
2010), we obtain conservation of volume of ocean water, i.e. 

Qin + Qout = 0. (3-4) 

With these two relations in Equations 3.3 and 3.4, we now consider the the “bulk 
mixing” model of river and ocean waters, adding freshwater discharge from a given 
river with salinity concentration Sr at river discharge rate Qr into the into the 
considered finite-volume cell, which gives an additional salinity flux into the box of 
SrQr, as illustrated in Figure 3-l(b). Next, we re-apply the conservation of water 
(volume) and of salt to this new situation. 

Starting first with conservation of volume (water mass), since the water mass flow 
rate into and out of the cell must still be equal (the incompressibility and neglected 
free-surface principle utilized to obtain Equation 3.4 still apply), we have Qr of 
additional water entering and also leaving the finite-volume V. 

However, considering the salt conservation, we also have a new salinity flux out 
of the cell, denoted SQr. This S is defined as the average salinity concentration 
of the entire finite-volume V and is different from Sin and Sout- This S represents 
the fully mixed salinity concentration within the cell. If the control volume becomes 
infinitesimally small, then the salinity concentrations become equal such that Sin = 
Sout = S. Using the average salinity concentration S models the exit of well-mixed 
salinity in the control volume. 

We now further define the mass of salt in the numerical finite-volume V as Ms = 
S -V. We can then simplify the left-hand side for a constant numerical finite-volume 
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I 

I ^outOput 

(a) 

Figure 3-1: Diagram of salinity fluxes into 
model without freshwater discharge and (b) 



(b) 


and out of the control volnme for (a) a 
i model with freshwater discharge added. 


V: 




and expand Equation 3.1 to include all four salt fluxes, to get 


dS • • 

^ ~jT SifiQifi + SjiQji SoutQout SQji. 

dt 


(3.5) 


(3.6) 


Now we can group terms with the same flow rate, recalling Equation 3.4 


dS 

dt 


IQ Q\ \ IQ Q \ Cottt 

yuR o)^ yOiji ^out) ^ 


(3.7) 


The first term in that final equation [Sr — S)^^ provides the additional rate of change 
in salinity that locally forces each finite-volume that is part of the coastal volume over 
which the bulk river mixing is applied. The left hand-side ^ represents the time- 
rate-of-change processes that affect salinity in that finite-volume and the last term 
[Sin — Sout)^^ represents the sum of all in and out fluxes that were here modeled as 
bulk terms. When the control volnme becomes infinitesimal, this last term becomes 
the snm of advective and diffnsive salinity flnxes; the first term, becomes a partial 
derivative with respect to time; and the additional river forcing term remains as is, 
i.e. {Sr — S)^ excepted that ^ becomes the local mass flnx of river water (at a 
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point in the coastal ocean). 


From the above Equation 3.7, one then needs to distribute a priori where the river 
water is discharged in the coastal volume over which the bulk river mixing is applied; 
essentially, how to assign ^ in space. To do so, if we assume, for now, that each 
hnite-volume in that coastal volume receives the same amount of river water (the river 
water is equally distributed to each hnite-volume), and that all hnite-volumes have 
the same volume, then that total coastal volume over which the bulk river mixing 
is applied is Volume = NV where N is the number of hnite-volume cells used (per 
river), and V is still the volume of each hnite-volume cell. The mass of salt in that 
total volume is then Ms = S ■ Volume. In that case, the above equation becomes 

— = {Sr-S)^ + {Sin - Sout)^ (3.8) 


where Qr is now the total river discharge and Qout is the total advection hux. Of 
course, each hnite-volume cell does not have the same volume. In our actual discrete 
model, the river water is distributed according to the horizontal footprint of these 
hnite-volumes, which are all the same in the present MSEAS conhguration (Haley 
and Lermusiaux, 2010). 

If we ignore, for now, this "surface-based distribution of river waters" and work 
with the global Equation 3.8, we can utilize the Mirabito et ah (2012) dehnition for 
r: 


T = 


NV 

“ 5 " 


(3.9) 


we can rewrite the total fresh-water river forcing {Sr — S)^ (to be applied locally 
on each hnite-volume if that forcing is equally distributed in horizontal space) as 


{Sr - S) 
Tr 


(3.10) 


We now expand the salinity forcing function above to account for all the rivers in 
our model. Then the salinity forcing functions are inserted into the MSEAS primitive- 
equation ocean model. Specihcally, they are applied as a second term on the right 
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hand side of Equation 2.5: 


DS 


^River { O 

ai{x,y,z )^-^— 
i=l ^R,i 


S) 


(3.11) 


where Njnyer = H in the current iteration of our model, and ai{x, y,z) = 1 dehnes the 
coastal volume over which the bulk mixing model for river i is applied (a*(a;, y,z) = 0 
outside this volume). The time constants are provided as inputs to the MSEAS 
primitive-equation model. However, as mentioned above, the total river discharge is 
distributed based on the horizontal footprint of ai{x, y, z), and the vertical dependence 
will be discussed in Section 3.2.2. 


3.1.2 A Note on Indexing 

In Mirabito et ah, the index for tau is (j, k). This corresponds to river k in basin j 
and each unique {j, k) corresponds to a unique index i such that each river is assigned 
a unique index i. In Section 3.1.1, the generic subscript R is used to indicate that 
the terms are specific to any selected river and should be adjusted for each river i. In 
addition, Mirabito’s notation uses Q for discharge flow rate with units of m^/s. This 
thesis uses the symbol Q to indicate the same quantity using Q to denote discharge 
(or precipitation) amounts, as seen in Section 2.5. 


3.2 Potential Improvements to the River Discharge 
Model and Bulk Mixing Model 

The original river discharge model (Section 2.5) and bulk mixing model of river and 
ocean waters (Section 3.1.1) indicated improved forecasting of ocean parameters over 
the ocean model without the salinity relaxation, based on comparisons with data 
collected during the QPE experiment. Potential areas for future improvement were 
also identifled: time series estimation, depth of salinity relaxation, and number and 
shape of selected grid (finite-volume) cells were specifically sources of potential error 
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and future improvement (Mirabito et al., 2012). The following section reports research 
completed to address some of these areas. 

3.2.1 Time Series Estimation for River Discharge 

The river discharge model uses actual river discharge data for two of Taiwan’s ma¬ 
jor rivers, the Zhuoshuh and the Gaoping. The data was collected at hydrological 
stations near the mouth of each river that remained active during the period of 06 
-10 August 2009, when Typhoon Morakot passed over the island (Mirabito et ah, 
2012). The discharge rates of the two rivers are averaged together, then normalized 
and distributed among the major rivers in each basin using relative annual mean 
discharge values. See Section 2.5.2 for a more thorough description. Therefore, the 
river discharge model assumes that all rivers will reach their peak discharge at the 
same time. The only exceptions are the Zhuoshuii and the Gaoping Rivers; the time 
series for these rivers use their actual discharge data in order to improve accuracy. 

One possible method of improving the accuracy of this river discharge model is to 
use a different time series for each basin or even each river, based on the geographical 
position of the river and the path of the typhoon. In other words, the discharge time 
series would peak earlier for rivers on the east side of the island, which would have 
been the first to encounter the torrential rains of Morakot (see Figures 2-10 (a) and 
3-2). 

Using a best track position graphic from Joint Typhoon Warning Center, shown 
in Figure 3-2, we know that the eye of the storm made landfall just before 1800Z on 
07 August 2009, and the eye moved offshore just before 0600Z on 08 August 2009 
(Cooper and Falvey, 2009, p. 23). The storm center took approximately twelve hours 
to pass over the width of Taiwan. It can be assumed that the time difference between 
the peak rainfall rates on the eastern and western sides of the island is approximately 
twelve hours, and therefore the time difference between the maximum river discharge 
rate on the two sides could be on the order of twelve hours. 

The raw discharge rate data for the Zhuoshuii and the Gaoping Rivers is shown 
in Figure 3-3. From this data, it is shown that the time difference between peak 
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Figure 3-2: Historical Best Track for Typhoon Morakot. (Cooper and Falvey, 2009, 
p. 23) 


discharge rates of the two rivers is on the order of hours. The Zhuoshuii River is 
located on the west-central coast of Taiwan, and in drainage basin j = 2; the Gao- 
ping River is located on the southwest coast and in drainage basin j = 1. See Figure 
2-10(a). The distance between the two rivers is approximately 154 km, or 96 miles 
(National Hurricane Center, 2014). Importantly, the duration of the large discharge 
is of the order of four days for both rivers. This duration is much larger than the few 
hours offset of the two peak river discharges. 

For our purposes, the time difference between peak discharges is also relatively 
small when compared to the timeframe of our ocean primitive-equation model, which 
forecasts for days and even weeks. Within the context of these mid-range forecasts, 
twelve hours’ difference is likely not sufficient to have an appreciable improvement in 
model accuracy. If the desired forecast period was very short-range, then the increase 
in model complexity may be compensated by the slight increase in accuracy. For this 
work, it was determined that forecast improvements from varying time series would 
be negligible. 

Improvements to the time series accuracy could be made by incorporating addi- 
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Figure 3-3: Discharge Data for the Zhuoshuii and Gaoping Rivers from 5-14 August. 
(Mirabito et ah, 2012) 

tional discharge data into the river discharge model. This is discussed in Chapter 

5. 

3.2.2 Depth of Selected Finite-Volume Grid Cells in the Bulk 
Mixing Model 

In the original bulk mixing model of river and ocean water, the salinity was relaxed 
to a uniform depth of 5 m (Mirabito et ah, 2012). This depth was chosen to represent 
the depth of the river discharge. Since the relaxation salinity models a well-mixed 
average salinity within a selected grid cell, there is no constraint on applying the 
salinity relaxation only near the surface of the selected grid cells. Applying the 
salinity to a lower depth would model vertical mixing. During the passage of a large 
typhoon such as Morakot, it can be expected that vertical mixing in the upper layers 
of the ocean would occur. Subsequent improvements to the original model used a 
depth of 10 m, which was the starting point for this study. 

In order to determine if 10 m is the appropriate depth for the salinity relaxation, 
we studied the data collected soon after the passage of the typhoon by the QPE 
experiment. 

The OR2 data, shown in Figure 3-4, was taken off the northeast coast and near 
the mouth of the Lanyang River (index i = 8). We look at coastal station 39, as well 
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(b) 

Figure 3-4: (a) Map of CTD drop locations northeast of Taiwan, taken by the R/V 
Ocean Researcher 2 (OR2) during the period 13 - 17 August, 2009. (b) Salinity 
profiles from these CTD casts. (Leslie, 2012) 

as stations 33, 34 and 38. To see indications of freshwater transport away from the 
coast, we use stations 1, 31, 32, 35 and 40. At station 39, we see low salinity, less 
than 33.8 psu, to a depth of approximately 50 m. The low salinity, less than 33.6 
psu, is present to a depth of approximately 20 m. Stations 32 - 36 also have very low 
salinities, less than 33.6 psu, to a depth of approximately 10 - 15 m. 

In the OR3 data, shown in Fignre 3-5(b), we have data jnst offshore for rivers 
with indices 3 — 6. Stations 1 - 3, 7 - 10, 23 - 25, and 34 - 37 are all close to shore 
and could reflect lower salinities from freshwater discharge. In addition, stations 12, 
26, 33 and 38 may indicate the transport of this freshwater from near-shore waters. 
Upon examination of the data, we see stations 1-3 have very low (less than 33.6 psu) 
salinities to a depth of approximately 15 m. Stations 7-10 show lowest values to the 
ocean bottom, which is approximately 15 m. Stations 23 - 25 show salinities below 
33.6 psu to a depth of about 25 m, and stations 34 - 37 indicate lowest salinities only 
to about 15 m. 

In general, the salinities are much lower on the western side of the island than on 
the northeast side, typically below 33.9 psu while the northeast casts show salinities 
above 34.4 psu. Per Newhall et al. (2010), the water mass northeast of Taiwan has 
typical salinity ranges of 33 - 34.8 psu. After the passage of Typhoon Morakot, 
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Figure 3-5: (a) Map of CTD drop locations northeast of Taiwan, taken by the R/V 
Ocean Researcher 3 (OR3) during the period 13 - 17 August, 2009. (b) Salinity 
profiles from these CTD casts. (Leslie, 2012) 

“typhoon-induced brackish warm water was found off the northen tip of Taiwan. The 
salinity was as low as 32 around the northwestern tip of Taiwan. Since the runoff 
was modest for the rivers in the northern Taiwan, this freshwater souree must have 
originated from the southwest coast of Taiwan.” The oceans west of Taiwan are 
over the shelf, shallower and fresher than the slope waters to the northeast, where 
salty Kuroshio water mixes with the shelf waters flowing down the slope. Specific 
profiles were identified by Newhall et al. (2010) as being impacted by freshwater river 
discharge. These salinity profiles and their locations are shown in Fignre 3-6. 

From the above analysis, several new test depths were selected. Depths of 20, 35 
and 50 m were selected to test the extremes and to allow a bracketing of the best 
mixing depth. Recall from Eqnation 3.9 that the relaxation time-scale T{t), where t 
is time, in the salinity forcing Eqnation 3.10 depends on the volume of the selected 
finite-volume cells, which is determined by the number of cells and the depth to which 
the salinity relaxation factor is applied, i.e. 

NV 

T = —^ 

Q 

The time series of r depends on the volume over which the salinity relaxation is 
applied, and was therefore re-compnted for each of the new test depths 20, 35 and 50 
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T(z) and S(z) are 
plotted for the 
stations outlined in 
blue. S < 33 is 
found in green 
regions. S < 32.5 is 
found in red area. 
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Figure 3-6: Salinity profiles and locations of low salinity observations in the wake of 
Typhoon Morakot. (Leslie, 2012) 
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m and applied to the existing river model. After these initial tests and our analysis of 
their results, as well as updates to the area over which we apply the bulk mixing, i.e. 
the hnite-volume grid cell conhguration, and changes to the grid cell conhguration 
discussed in Section 3.2.4, new model runs to depth of 10, 15, 20 and 35 m were 
completed. See Chapter 4 for a more detailed explanation. 

3.2.3 Salinity Relaxation Values 

The original river discharge model tested multiple values of river salinity Sr; 10, 15 
and 20 psu were tested. The best choice was selected by comparison with observations, 
and Sr = 10 psu was deemed the most accurate (Mirabito et ah, 2012). Subsequently, 
data was obtained to indicate that Sr = IS psu gave more accurate results when 
compared with QPE data, and with other minor changes to the model. For this 
study, Sr = 18 psu was used. 

If the model output salinity approaches the river, or relaxation, salinity (S and Sr, 
respectively, in Equation 3.10), then the river discharge model is no longer approxi¬ 
mating an input of relatively fresh water as compared to the surrounding seawater. 
The model output is confirmed to maintain minimum salinities well above the relax¬ 
ation salinity of 18 psu. 

3.2.4 Number and Shape of Selected Grid Cells 

The number and configuration of grid cells for salinity relaxation was another possi¬ 
ble area for significant improvement of the bulk mixing model. The original model 
was tested with n=4 and n=7 grid cells per river (Mirabito et ah, 2012). However, 
simulations with this bulk model indicate a too small foot print for the rivers and 
too low salinities when compared to observed values near the river mouths. Hence, 
hindcast simulations using a larger footprint were determined to more closely match 
the available data. 

Hence, six new configurations, or shapes, were created and tested. All used the 
original n=7 grid cells per river and therefore used the original time series. The shape 
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of the configurations was based off surface salinity model forecasts from the original 
river discharge model; the grid cells selected for salinity relaxation were chosen to 
mirror the shape of the low salinity areas near the mouths of the rivers. A sampling 
of new configurations is shown in Figure 3-7. 

After evaluating these results, it was found that even n=7 may not be a large 
enough footprint to accurately capture the bulk mixing occuring from the freshwater 
discharge input. Three more configurations, consisting of n=10, n=14, and n=25 
finite-volume grid cells per river mouth were created and tested. The large number of 
cells and the close proximity of the river mouths, especially for rivers z = 4 — 6 (the 
Wu, Dajia and Daan Rivers, respectively) resulted in one large combined footprint 
on the western side of the island. These configurations are shown in Figure 3-7. 

3.3 Sensitivity Studies of the MSEAS Ocean Model 

Section 3.2 discusses many updates that could potentially improve the accuracy of 
the river discharge model and bulk mixing model. Since the river discharge model is 
used as an input into the MSEAS primitive-equation ocean model, it would be remiss 
to ignore sensitivity studies of this larger model. One such study done during the 
scope of this project was a "restart" study. This method was used to examine the 
effects of tidal forcings and atmospheric forcings on oceanographic properties. 

The concept is to run the model for a pre-determined time, then remove a forcing 
and let the model continue to run. A “control” run is also done, keeping the forcings 
in place throughout the run. Then, a comparison can be made between the restart 
and the control to determine the effects of the forcing. 

This study included an analysis of the impacts of the atmospheric forcings and 
the tidal forcings. Two sets of restarts were run, with each one removing one of the 
forcings. The forcings were removed after 7, 14, 21 and 28 days; the total time- 
frame for this study was 37 days. Each restart run was completed with and without 
data assimilation; the assimilation scheme was adapted to compensate for the new 
timeframes. 


54 






111 





1 


|pn 


,■ ip 


nir 

IvFTV 

ii: 


(d) (e) (f) 


Figure 3-7: Comparison of grid cell configurations. Salinity relaxation grid cells are 
marked with red boxes, (a) Original configuration, (b) A new configuration, called 
Test 5, maintaining original n=7 grid cells per river, (c) Another new configuration, 
called Test 6, with n=7 grid cells per river, (d) A configuration, called Test 7, 
increasing to n=10 grid cells per river, (e) A configuration, called Test 8, increasing 
to n=14 grid cells per river, (f) A configuration, called Test 9, increasing to n=25 
grid cells per river. 
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Chapter 4 


Results 


We study the results of the different runs to determine the best combination of pa¬ 
rameters for accurate salinity forecasting. All Figures in this section show model runs 
with no data assimilation, unless otherwise noted. 


4.1 Effects of Footprint Shape on the Bulk Mixing 
Model 

We begin by examining the results of MSEAS primitive-equation ocean model runs 
with no river forcing and with the original river forcing from Mirabito et al. (2012). 
We then compare these to new runs done with different overall shapes in the grid 
cell footprint. All of these first runs did not utilize any changes in the total volume 
over which the salinity relaxation was applied; that is, the volume was n=7 and 10 
m depth for these initial tests. 

We note that a large number of simulations were completed. Here, we only report 
on a subset of these simulations, selected to illustrate the results. 

Comparing the surface salinity forecasts from the the MSEAS model run without 
river forcing (Eigure 4-3) to the salinity forecasts utilizing the original river discharge 
model (Eigure 4-4) shows the marked difference made by the use of the river discharge 
model. The lower salinities in the coastal waters off Taiwan, especially on the western 
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QPE OR2 and OR3 Data 



Longitude 

Figure 4-1: Map of CTD drop locations taken by the R/V Ocean Researcher 2 (OR2) 
dnring the period 13 - 17 Augnst, 2009. This is a combination of the data prohles 
shown in Figures 3-5(a) and 3-4(a). (Leslie, 2012) 

side of the island, more closely match the data obtained dnring the QPE DRI. Fignres 
4-5 and 4-6 show the surface salinity forecast for our Tests 5 and 6, respectively. Each 
keeps the original n=7 grid cells, but in a different geographical configuration than 
the original model (see Eigure 3-7(a)-(c)). The different shapes were created with 
the overall ocean flow in mind; the shapes were drawn to show strong advection 
into surrounding currents, spread along the coast, or discharge directly away from 
the coast, ignoring ocean flows. Despite the different shapes, the snrface salinity 
forecasts appear similar to the original river model. 

In Eignre 4-7(a), (b), and (c), we show the differences between these forecasts 
and the data. Any differences are very snbtle. These plots inclnde a profile for all 
available data, over 100 profiles, corresponding to CTD cast locations shown in Figure 
4-1. Our data is for 13 - 17 August 2009, which is only a few days after the passage of 
Typhoon Morakot. Thus, we do not expect to see significant salinity impacts at large 
distances from the coast. Having too much irrelevant data can disguise our results; 
hence, we choose to focus only on the profiles that would be affected by freshwater 
river discharge. In particnlar, six profiles have been identified as being impacted by 
river discharge (Newhall et ah, 2010). We now generate plots comparing the model 
forecast salinities to the salinity data for only these plots. Using this approach, subtle 
but significant differences between model runs become more obvions. 
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QPE River Influenced Profiles 



QPE River Influenced Profiles 
identifiers 



(a) 


(b) 


Figure 4-2: Salinity profiles taken 13 - 17 August 2009 for the six locations that have 
been identified as being impacted by freshwater river discharge by Newhall et al. 
(2010). (a) shows salinity variations with depth in these six profiles, while (b) shows 
the locations of these profiles. 


We now study Figure 4-8, which compares the river-impacted profiles for the 
original model, as well as our new Tests 5 and 6. This Figure shows that the impact 
on salinity in river-impacted profiles is still very similar for the three forecasts. There 
is very little difference evident between plots of different cell configurations when the 
configurations have the same number of cells; however, Test 6 is slightly better than 
the other two options. 


4.2 Effects of the Horizontal Extent of the Bulk Mix¬ 
ing Model 

We design two new configurations, increasing the the size of the fresh-water footprint 
from a number of cells n=7 to n=10 and n=14. We also generated a very large 
footprint model. Test 9, with n=25 grid cells. These shapes are diagrammed in 
Figure 3-7(d)-(f)). As before, the various shapes were carefully created to account 
for ocean currents, studied from runs without salinity forcing, as well as to mirror 
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salinity mixing patterns seen in previous runs. We chose not to investigate a smaller 
footprint, as this was previously studied by Mirabito et al. (2012) and determined 
less accurate than the original n=7 configuration. 

The comparisons of these tests are shown in Figure 4-7(d)-(f). Again, the full 
profile plots obscure the desired data, so we analyze the plots with selected river- 
impacted profiles in Figure 4-9. 

Here, we see a clear improvement using the n=25 configuration. There is a reduc¬ 
tion in salinity differences from the data in all profiles. 


4.3 Effects of the Depth of the Bulk Mixing Model 

A separate study was done to determine the effect of the salinity relaxation depth on 
the model salinity forecasts. The original model, as well as the new configurations, 
were run with different bulk mixing model depths and the results compared in Figure 
4-10. Again, the differences from the data were plotted for all the profiles, as well as 
for the river-impacted profiles. 

As seen in the the previous results, Test 6 was slightly better than the original and 
Test 5 configurations (Figure 4-11). Comparing Test 6 with the forecasts for larger 
footprints (Figure 4-12) again shows that Test 9, with n=25 cells, is slightly better 
than the other tests, although the improvement is less marked at 20 m than at 10 
m. The configurations were also run using a salinity relaxation depth of 35 m, which 
showed less improvement than the 20 m and is not included here. 

As the n=25 configuration was the best choice at both 10 and 20 m, this new 
configuration was run with another, new salinity relaxation depth of 15 m. Then 
these three runs were compared to the data; the comparison of results is shown in 
Figure 4-13. The combination of the largest footprint and the original depth most 
closely matches the data from 13-17 August. 
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4.4 Final Results 


Figure 4-14 shows the surface salinity hindcast for the MSEAS primitive-equation 
ocean model run using the n=25 configuration with salinity relaxation to a depth 
of 10 m. Figure 4-15 shows the same model run, using data assimilation techniques 
(Haley and Lermusiaux, 2010) (Haley et ah, 2014). 

On August 7, the current pattern shifts on both the western and eastern side of 
the island, following the northwestward track of Typhoon Morakot (refer to Figure 
3-2). Beginning on August 9, the effects of the river discharge become obvious with 
the low salinities concentrated in the coastal waters around Taiwan. This is most 
noticeable on the western shore, where 4 large rivers flow into the sea relatively close 
together. These salinity lows are evident throughout the forecast period, although 
they are strongest on August 9, in the middle of the period of heaviest rainfall. 
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(a) August 07 


(b) August 09 
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(f) August 25 


Figure 4-3: Surface salinities forecast by the MSEAS ocean model with 
forcing. This run was completed for comparison purposes. 
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(c) August 13 


(d) August 17 



Figure 4-4: Surface salinities forecast by the MSEAS ocean model with the original 
river forcing, using n=7 finite-volume cells, with salinity relaxation applied to a depth 
of 10 m. The configuration is shown in Figure 3-7(a). This run was completed for 
comparison purposes. 
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Figure 4-5: Surface salinities forecast by the MSEAS ocean model with configuration 
Test 5, using n=7 finite-volume cells, with salinity relaxation applied to a depth of 
10 m. The configuration is shown in Figure 3-7(b). 
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(e) August 19 (f) August 25 

Figure 4-6: Surface salinities forecast by the MSEAS ocean model with configuration 
Test 6, using n=7 finite-volume cells, with salinity relaxation applied to a depth of 
10 m. The configuration is shown in Figure 3-7(c). 
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(a) Original (7 cells) @ 10 m 



(c) Test 6 (7 cells) @ 10 m 







(e) Test 8 (14 cells) @ 10 m 



(b) Test 5 (7 cells) @ 10 m 



(d) Test 7 (10 cells) @ 10 m 



(f) Test 9 (25 cells) © 10 m 


Figure 4-7: Comparisons of salinity data with MSEAS salinity forecasts using different 
configurations of the bulk mixing model, with salinity relaxation to a depth of 10 m. 
Plots show the data subtracted from the specified model output. The profile locations 
are shown in Figure 4-1. 
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(a) Original (7 cells) @ 10 m 



(b) Test 5 (7 cells) @ 10 m 




(c) Test 6 (7 cells) @ 10 m 

Figure 4-8: Comparisons of salinity data (in profiles with known river discharge in¬ 
fluence (Newhall et ah, 2010)) with MSEAS salinity forecasts using different configu¬ 
rations of the bulk mixing model, with salinity relaxation to a depth of 10 m. Plots 
show the data subtracted from the specified model output. The profile locations and 
data are shown in Figure 4-2. 





















(a) Test 6 (7 cells) @ 10 m 



(b) Test 7 (10 cells) @ 10 m 



(c) Test 8 (14 cells) @ 10 m 


(d) Test 9 (25 cells) @ 10 m 


Figure 4-9: Comparisons of salinity data (in profiles with known river discharge in¬ 
fluence (Newhall et ah, 2010)) with MSEAS salinity forecasts using different configu¬ 
rations of the bulk mixing model, with salinity relaxation to a depth of 10 m. Plots 
show the data subtracted from the specified model output.The profile locations and 
data are shown in Figure 4-2. 
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(c) Test 6 (7 cells) @ 20 m 


(d) Test 7 (10 cells) @ 20 m 



(e) Test 8 (14 cells) @ 20 m 


(f) Test 9 (25 cells) © 20 m 


Figure 4-10: Comparisons of salinity data with MSEAS salinity forecasts using dif¬ 
ferent configurations of the bulk mixing model, with salinity relaxation to a depth of 
20 m. Plots show the data subtracted from the specified model output. The profile 
locations are shown in Figure 4-1. 
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(a) Original (7 cells) @ 20 m 
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(b) Test 5 (7 cells) @ 20 m 




(c) Test 6 (7 cells) @ 20 m 

Figure 4-11: Comparisons of salinity data (in profiles with known river discharge 
infiuence (Newhall et ah, 2010)) with MSEAS salinity forecasts using different config¬ 
urations of the bulk mixing model, with salinity relaxation to a depth of 20 m. Plots 
show the data subtracted from the specified model output. The profile locations and 
data are shown in Figure 4-2. 
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(a) Test 6 (7 cells) @ 20 m 
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(b) Test 7 (10 cells) @ 20 m 



(c) Test 8 (14 cells) @ 20 m 


(d) Test 9 (25 cells) @ 20 m 


Figure 4-12: Comparisons of salinity data (in profiles with known river discharge 
infiuence (Newhall et ah, 2010)) with MSEAS salinity forecasts using different config¬ 
urations of the bulk mixing model, with salinity relaxation to a depth of 20 m. Plots 
show the data subtracted from the specified model output. The profile locations and 
data are shown in Figure 4-2. 


71 




























Sa linity' Ortlerence: ModelOutput-alldata 13.7augrK/ 



Profile Nunrtier 


(a) Test 9 (25 cells) @ 10 m 
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(b) Test 9 (25 cells) @ 15 m 


Sa lin ity Ditterence: ModelOulput-all data 13^7augr Kr 



Profile Nunrfcer 


(c) Test 9 (25 cells) © 20 m 

Figure 4-13: Comparisons of salinity data (in profiles with known river discharge 
infiuence (Newhall et ah, 2010)) with MSEAS salinity forecasts using the same 25 
cell configuration of the bulk mixing model, with salinity relaxation to a depth of 10, 
15, and 20 m. Plots show the data subtracted from the specified model output. The 
profile locations and data are shown in Figure 4-2. 
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(e) August 19 


(f) August 25 
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Figure 4-14: Surface salinities forecast by the MSEAS ocean model with configuration 
Test 9, using n=25 finite-volume cells, with salinity relaxation applied to a depth of 
10 m. The configuration is shown in Figure 3-7(f). 
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(c) August 13 


(d) August 17 



(e) August 19 (f) August 25 

Figure 4-15: Surface salinities forecast by the MSEAS ocean model with configuration 
Test 9, using n=25 finite-volume cells, with salinity relaxation applied to a depth of 
10 m. These results include data assimilation. The configuration is shown in Figure 

3-7(f). 
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Chapter 5 


Conclusion 


This work focused on the improvement of ocean forecasting techniques in the highly 
dynamic and complex region surrounding Taiwan. We did this by focusing on a pe¬ 
riod of time immediately following the passage of Typhoon Morakot, which brought 
an unprecedented amount of rainfall within a very short time period. This freshwa¬ 
ter resulted in excessive river discharge, and the work of this thesis was to model 
this discharge via a bulk mixing and river discharge model to improve the ability of 
primitive-equation ocean models in forecasting salinity responses to extreme weather 
events. 

By varying the shape and size of the bulk mixing model footprint, as well as 
the depth, and examining the resulting impacts on ocean salinity forecasts, we were 
able to determine the optimal combination of salinity relaxation factors for highest 
accuracy. Utilizing data collected by the Quantifying, Predicting, and Exploiting 
Uncertainty Department Research Initiative (QPE DRI), the bulk mixing and river 
discharge model was carefully tuned to optimal performance. 

Understanding the impact of extreme weather events on oceanic properties and 
improving ocean model forecasting allows users to predict the changing environment 
more accurately. Enhanced knowledge of the environment gives the advantage to the 
decision maker: knowledge is power, and in an environment as vast and complex as 
the ocean, even a small advantage can make a big difference. 
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5.1 Future Work 


Research begets research, and this work is no exception. During the course of this 
project, many other areas of potential improvement were discovered. Some of these 
were considered for this work, but constraints of time prohibited a thorough inves¬ 
tigation of all possibilites. We include some of these here for others interested in 
furthering this work. 

5.1.1 Variation of Salinity Relaxation Factors 

The river discharge model and bulk mixing model, in its current form, applies salinity 
uniformly throughout a set number of finite-volume cells assigned to an individual 
river per Equation 3.11. In reality, we expect that the impact of the freshwater 
discharge is felt most strongly in ocean waters closest to the river mouth, with effects 
dissipating as distance from the river mouth increases. Factors such as annual mean 
discharge, topography, instantaneous river discharge rates, and riverbed depth create 
much variation between the individual river impacts on the surrounding ocean waters. 
Modeling these smaller effects may have significant impact on the river discharge and 
bulk mixing models’ ability to improve ocean forecasts. A trade-off between model 
complexity and accuracy must be considered. 

The number of finite-volume grid cells per river could be adjusted to a variable 
value, with each river having a different size footprint within the model. The salin¬ 
ity relaxation could also be applied at variable depths within each river’s footprint. 
The salinity value used in the bulk mixing model could be applied non-uniformly 
throughout the footprint. 

5.1.2 MSEAS Primitive-Equation Model 

This study utilized the MSEAS primitive-equation ocean model. The many different 
model runs were run over the same domain as standalone runs. A study of the 
impact of domain size and/or grid resolution would enhance our knowledge of the 
model’s response to the river discharge model. The model is capable of utilizing 


76 



nested runs to increase resolution, and this technique would also contribute to a 
better understanding of the interactions between the primitive-equation model and 
the river discharge model. 

5.1.3 New Discharge Data 

Since the development of the original river discharge model, more river discharge 
data have become available. Jan et ah (2013) has discharge data for two additional 
rivers, the Zengwen (i = 2) and the Wul (i = 4) Rivers (Figure 5-1). This data 
could help to develop a more accurate time series to be used for the river discharge 
model. An interesting feature of this data is that the Zengwen River has a peak 
discharge rate that is approximately equal to the peak discharge rate of the Gaopmg 
River. Recall from Table 2.1 that the Zengwen River has an annual mean discharge 
of 50 m^/s, while the Gaopmg River has an annual mean discharge rate of 220 m^/s. 
The Wu River, with an annual mean discharge of 120 m^/s, has a peak discharge 
rate that is an order of magnitude lower than that of the other three rivers. This 
indicates that dividing the discharge between rivers based on annual mean discharge 
rates introduces a source of error. However, Figure 2-10(a) shows the geographical 
location of these four rivers. The Wu River is the furthest north; the other rivers 
are sourced high in the mountains of the southern part of Taiwan. It is important to 
note that the heaviest rainfall ocurred on the southern part of the island (see Figure 
2-8). Then it would be expected that the Gaopmg, Zengwen, Zhuoshuii Rivers would 
have higher peak discharge rates from Typhoon Morakot rainfall; this supports the 
basin-averaged rainfall distribution approach of Mirabito et ah (2012). 
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Figure 5-1: Discharge Data for the Gaoping (z = 1), Zengwen (i = 2), Zhuoshuii 
(z = 3), and Wu (z = 4) Rivers from August 2009. (Jan et ah, 2013) 
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